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I. INTRODUCTION 

The understanding of quark fragmentation functions has been rapidly evolving in recent 
years p]. Based on precise data on inclusive hadron production in hard scattering processes [2rl9], 
empirical fragmentation functions have been extracted, and their behavior under scale evolutions 
has been investigated in detail |10H12| . It is reasonable to expect that in the near future an 
understanding of the fragmentation functions will be attained which is as precise as the present 
knowledge on quark distribution functions. Very interesting prospects for the study of fragmenta- 
tion processes were opened by the HERMES [13] and JLab |14H17j semi-inclusive measurements 
on nuclear targets, which will be followed up by the planned Electron-Ion Collider |18} 119]. 

These exciting developments present a challenge for effective theories of QCD. Because both the 
distribution and fragmentation functions are basically nonperturbative quantities, effective quark 
theories are powerful tools for theoretical studies. In particular, the Nambu-Jona-Lasinio (NJL) 
model [20|, [21] , which was very successful in the description of quark distribution functions |22] , 
has been combined recently with the ideas of the jet model of Field and Feynman [23] to give a 
realistic description of quark fragmentation functions [23]. In recent work we have applied this 
NJL-jet model to the fragmentation functions for pions and kaons, and compared the results to the 
empirical functions [25]. The main advantage of this model, which is based on the multiplicative 
ansatz (product ansatz) of Field and Feynman, is that the fragmentation functions automatically 
satisfy the important momentum and isospin sum rules without introducing any new parameters 
into the theory. Ultimately, the NJL-jet model will provide a consistent framework to describe 
semi-inclusive deep inelastic lepton-hadron scattering, as well as allowing predictions for nuclear 
targets. 

The purpose of the present paper is to extend the NJL-jet calculations of quark fragmentation 
functions in the following three directions: First, we include the effects of secondary pions and 
kaons, which come from the decay of intermediate p, K* , and <f> mesons. We also include the 
strong decays of the vector mesons to 7r and K two-particle final states. Our aim is to investigate, 
in particular, the role of the p meson for the softening of the pion momentum distribution. The 
effects of three-particle decays, like the oj meson, will be left for future work, as those processes 
require a treatment of nontrivial phase space factors that cannot be evaluated analytically, which 
goes beyond the usual convolution formalism. 

Second, we will include the fragmentation processes to nucleons and antinucleons. For this pur- 
pose, we will use the splitting functions obtained in the quark-diquark description of baryons [26] , 
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taking into account only the effects of the scalar diquark as a first step in extending the model. The 
need to include the axial-vector diquark to fully describe the nucleon structure and fragmentation 
functions has been demonstrated in the earlier work |27H29| . thus this remains a priority for the 
future developments of the model. The NJL description of nucleons as bound states of a quark 
and a scalar diquark has already been applied to fragmentation functions in previous work [30} [31] . 
In the present paper we will, however, go beyond those earlier attempts by including the elemen- 
tary fragmentations to nucleons and antinucleons in the cascadelike processes also including pions, 
kaons and vector mesons. 

Third, we use the Monte Carlo (MC) method to solve for the fragmentation functions within 
the quark-cascade model as opposed to solving integral equations in previous work. We demon- 
strate the viability of this approach to replace the integral equations by providing very similar 
order of precision in determining the fragmentation functions, at the same time allowing to relax 
the approximations necessary for formulating the integral equations and easily incorporating the 
resonance decays i nto the model. 

This paper is organized as follows: In Sec. [XTJ, we present the calculations of the vector meson 
elementary fragmentation functions and elementary fragmentation functions of a quark to nucleon 



antinucleon pair. In Sec. Ill, we describe the Monte Carlo method for calculating fragmentation 



functions and include the vector meson decays. In Sec. IV we present the resulting solutions for 
the fragmentation functions. Section |V| contains some concluding remarks and a future outlook for 
extending the model. 



II. ELEMENTARY FRAGMENTATION FUNCTIONS 



In this article, we extend our previous work of Ref. [25J using the same notation and model 
parameters. In the current section, we evaluate the "elementary" fragmentation functions of quarks 
to hadrons as a "one-step" process in the NJL model using light-cone (LC) coordinate^ The NJL 
model we use includes only four-point quark interaction in the Lagrangian, with up, down, and 
strange quarks and no additional free parameters (see, e.g. Refs. [32-34J for detailed reviews of 
the NJL model). We employ Lepage-Brodsky (LB) "invariant mass" cut-off regularization for loop 
integrals (see Refs. [53] for a detailed description as applied to the NJL-jet model), except when 
calculating meson-quark couplings as discussed further in this section, and use our previous values 
for the constituent quark masses for light and strange quarks M u = 300 MeV, M s = 537 MeV. 

1 We use the following LC convention for Lorentz 4- vectors (a + ,a~,ax), a = -jn{ a ° ± a 3 ) and ai = (a 1 , a 2 ). 
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A. Pseudoscalar Meson Splitting Functions 
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FIG. 1. Quark splitting function for mesons. 

The elementary fragmentation function of quark q emitting a meson m carrying light-cone 
momentum fraction z is depicted in Fig. [I] For completeness of the description we present the 
results for the pseudoscalar mesons derived in [25], leaving out the details. In the frame where 
the fragmenting quark has kj_ = 0, but nonzero transverse momentum component = —~pa_/z 
with respect to the direction of the produced hadron, the relation for the elementary fragmentation 
function is 

C m f d^k 

W = -^9 2 mqQ \ J j^Tr[S 1 (k) 7 + S 1 (k) l5 ^-p + M 2 ) 75 ] 

xS(k- - p-/z)2n5((p- kf - Mf) (1) 

_C7 2 [ d? P± ^ + (( z -i) Ml + M 2 ) 2 

2 9rnqQZ J (2vr) 3 (p2 + z(z - l)M? + zM$ + (1 - z)m 2 J 2 ' 1 ' 

Here Tr denotes the Dirac trace and the subscripts on the quark propagators denote quarks of 
different flavor - also indicated by q and Q, where the meson of type m under consideration has 
the flavor structure m = qQ. The C™ is the corresponding flavor factor given in the Table |lj The 
pseudoscalar meson-quark coupling constant, g mq Q, is determined from the residue at the pole in 
the quark-antiquark t-matrix at the mass of the meson under consideration. 

In LB regularization the cut-off in the transverse momentum, P? 3 is given by 

Pi = *(1 " *) (V A i + m2 m + V A 3 + M l) -(l-z)™ 2 m -zMl (3) 

where A3 = 0.67 GeV denotes the 3-momentum cutoff obtained in Ref. [35], yielding the following 
values for the meson-quark couplings: 

9-nqQ = 3.15, g Kq Q = 3.39. (4) 

A consequence of LB regularization is a limited range for z in corresponding regularized func- 
tions, < z m i n < z < z max < 1, where z m i n and z max are determined by imposing the condition 
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P£ > in Eq. Q3p . These range limitations depend on the masses of hadrons and quarks involved. 
For example, the z limits are very close to endpoints (z = and z = 1) for splitting functions to 
pions, but are quite far from them for heavier hadrons like nucleons or eft meson. 



B. Vector Meson Splitting Functions 

In the experimental measurements of the quark hadronization process, one usually directly 
detects only the relatively long lived particles: pions, nucleons, and sometimes kaons. These 
particles can come either as direct emissions of the fragmenting partons or as decay products of 
hadronic resonances. We note that in the current article we will only consider the strong decays 
of the hadrons. The inclusion of the vector mesons in the NJL-jet model is important, since it 
has a two-fold effect on the description of the previously calculated pion and kaon fragmentation 
functions. First, the availability of the new emission channels decreases the normalization of the 
direct quark splitting functions to pions and kaons. In fact, in the early models for the quark 
fragmentation, it was assumed that the ratio of the vector to pseudoscalar meson fragmentations 
should follow from the simple spin state counting rule, 3:1. Later, it was shown experimentally that 
this vector meson ratio is significantly smaller, especially in the light quark sector (see, for example, 
Ref. [36j ) , and thus has been attributed to dynamic effects like the masses of the considered mesons 
having a strong influence on the ratio. (The above mentioned rule holds only for those mesons 
containing the heaviest quarks, where the differences in masses between pseudoscalar and vector 
mesons are negligible.) Second, the strong decays of the vector mesons produce pseudoscalar 
mesons with a z distribution that is distinguishable from the ones emitted directly in the quark 
fragmentation process, thus modifying the final z dependence of the fragmentation functions. 



The basic NJL interaction Lagrangian relevant for the vector meson channels has the form 
— G v Oip^X a i^) 2 , where a = 0,1,. ..8 with Ao = \/2/3 1. If the 4-Fermi coupling constant G v 
is chosen to reproduce the mass of the p meson as the pole of the qq t-matrix, the calculated 



masses of K* and 4> mesons agree well with the experimental values (see Appendix C 3 ) . Using the 
vector meson - quark vertex gmqQl^^a (versus tg m qQ7 5 Xa for pseudoscalar mesons), the elementary 
fragmentation function can be written as 



C 7 f d^h 

W= E -t^Qo / 7^4^[ 5 i( fc )7 + 5i(^(^-^ + M 2 ) 7 ,] e ^A,p)6-(A,p) 

A==Fl,0 J ^ ' 

x6{k- -p-/z)2Ti5({k-p) 2 -Ml), (5) 
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where the sum is over the polarization^] of the vector meson m and Tr is over Dirac indices. 
The details of the calculation of the Tr and integration over plus and minus components of the 
momentum are given in Appendix [B| Here we simply present the resulting expression: 

C™ 2 f d2 p±2 (pi+((l- 2 )M 1 -M 2 ) 2 )+^(( P i- Z 2 M 1 M 2+ (l- Z )m^) 2 +pi^(M 1 +M 2 ) 2 ) 



-9 mqQ Z J ^3 (^ +2 (^l)M2+^M| + (l- 2 )m^) 2 

(6) 

The vector meson-quark couplings, g m qQ, are determined from the residue at the pole in the 
quark-antiquark t-matrix at the mass of the meson under consideration. Here we encounter a 
deficiency in LB regularization scheme, where the mesons with mass larger than the sum of con- 
stituent quark's masses are not bound. This problem is usually circumvented [37J by choosing 
a large enough constituent quark mass accounting for the masses of all the considered hadrons. 
Here to keep consistency with our previous calculations, we choose to use a different approach: we 
will use the same constituent quark masses as in the previous NJL-jet model calculations, but in 
order to assess the vector meson-quark couplings we will use the proper-time (PT) regularization 
scheme of Refs. [371 - 139] that mimics confinement and eliminates the unphysical decay thresholds. 
The details of these calculations are presented in Appendix [Cj where we compare the pseudoscalar 
meson-quark couplings calculated both in LB and PT regularization schemes and ensure that they 
are practically the same. Also, the strange constituent quark mass determined from the experi- 
mentally measured kaon mass is practically the same in both regularization schemes. It is therefore 
reasonable to use the PT scheme for the calculation of the vector meson- quark couplings in the 
present model, while keeping the LB scheme for the regularization of the elementary fragmentation 
functions. 



C. N and N fragmentation channels 

In this section, we consider the elementary splitting functions, which involve nucleons and 
antinucleons, i.e., q — > ND and the subsequent process D —> Nq. Here q = u, d denotes the 
nonstrange quark, N = n, p denotes the nucleon, and D denotes a scalar diquark ( J = T = 0, 3 C ). 
The first process leads to the elementary fragmentation functions (z), represented as a cut- 
diagram in Fig. [2^,, and d®(z), which is shown in Fig. [2^. We recall that for a general process 
o — > be, the fragmentation function d h a {z) corresponds to the probability distribution of the LC 

2 Since in this paper we consider only spin-independent fragmentation functions, we include the spin degeneracy 
factor 2sj, + 1 into the definition of the fragmentation function d h a (z), so as to facilitate the comparison to the 
empirical parametrizations [101 112| . (We note that the operator definition used in Ref. |24| does not include this 
degeneracy factor.) 



(a) 



(b) 




(c) 



FIG. 2. Two contributions to the quark fragmentation function to a nucleon and antidiquark are depicted on 
Figs, (a) and (b), and the quark distribution function in nucleon fq(x) is depicted in (c) as a cut diagram 
(left) and as the equivalent Feynman diagram (right). 



momentum carried by b relative to a, and includes a sum over the spin-color quantum numbers of 
the spectator c. Therefore, the two fragmentation functions of Fig. [2] are simply related by 

df(z)= 1 -d^(l-z). (7) 

The second process D — > Nq leads to the elementary fragmentation functions d^(z), shown in 
Fig. j3^i, and (£L(z) of Fig. |3jo. The relation for this case is 

d^(z) = Jdg (1 - z) . (8) 

To evaluate the functions (z) and d^(z), one can either directly consider the cut diagrams of 
Figs. [2^i and[3^, or one can make use of the fact that all elementary NJL fragmentation functions are 
formally related to the more familiar distribution functions f(x) by crossing and charge conjugation, 
which is expressed by the relatiorj^] [3U1 HI] 

d b a (z) = (-l) 2 (-+^)+ 1 (2s b + I)- f b a (x = -), (9) 



la 



where j a is the spin-color degeneracy of a. Therefore, in order to obtain d^(z), we can use the 
well-known expression for the quark distribution function inside a nucleon, where the nucleon is 
described as a bound state of a quark (mass M) and a scalar diquarkj^] (mass Mn). This distribution 

3 We emphasize that this so called "Drell-Levy-Yan relation" can be used only to obtain the expressions for the 

elementary unregularized fragmentation functions from the distribution functions. 

4 The basic NJL interaction Lagrangian in the scalar diquark (D) channel has the form 

Gs {tpToCT 2 l3 A tp T ) (ip T C- 1 j 5 T 2 p A ip), where /3 A = ^3/2 \ A {A = 2,5,7) and C = 17270- The diquark 
mass Md is calculated as the pole of the qq t-matrix, and the nucleon mass as the pole of the qD t-matrix |39j . 
The 4-Fermi coupling constant G s is chosen so as to reproduce the experimental nucleon mass. 



s 



function fq{x) is shown in Fig. [2J3. It is easily evaluated by using the form of the nucleon vertex 
function 



Tn(p) = \ i-Z N —^-u N (p) 
V V- 



(10) 



where the Dirac spinors are normalized as itjvujv = 1, and the normalization factor Zjq is given by 



N 



(11) 



--M N 



with the quark-diquark bubble graph expressed by the Feynman propagators of the quark and the 
diquark: 

d 4 k 



/d 4 fr 
S F (k)A s (p-k), 

which yields for the normalization factor 



(12) 



dx 



dV 



(1 - x) (p 2 ± + (M N x + M) 2 ) 



J ( 27 [M%x{l — x) — p]_ — M 2 — x(M% - M 2 )] ■ 



(13) 



We obtain from the Feynman diagram of Fig. [2 

d 4 fe 



fq{x) = iM N Z N u N 



(2tt)' 



(S F (kh + S F (k)) A s {p - k)6(k^ - p-x)u N 



\z N {l-x) J 



d 2 k 



T 



k\ + (M N x + Mf 



l2 - 



(14) 



2 'J (2vr) 3 [ fc 2 + M 2 (l-x) + M 2 D x-M 2 N x{l-x)\ 

Here, (q,N) = (u,p) or (d,n), while for the other flavor combinations the distribution function 



vanishes. The relation (|9|) then gives 



d*(z) = ±Z N (l-z) f 



dV 



p\ + (M N + Mzf 



(15) 



(2tt) 3 [ p 2 _ M 2 z(1 _ z ) + M 2 Z + _ z )] 2 ' 

The function d®(z) can then be obtained from (j7|, which completes the evaluation of the fragmen 
tation functions of Figs. [2^ and [2)3. 



(a) 



(b) 



FIG. 3. Two contributions to the antidiquark fragmentation function to an antinucleon and a quark are 
depicted on Figs, (a) and (b). 
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For the fragmentation function d^(z) of Fig. 3i, we note that charge conjugation invariance 
gives d~(z) = dp(z), where the latter function is related to the diquark distribution function inside 
the nucleon (fp(x)) by the crossing relation (|9j). The function fp(x) in turn is simply obtained 



from the quark distribution function (14) by replacing x — > 1 — x. Summarizing, we obtain 

d»(z) = z~f q N (l-x)\ x=1/z 

-\ Z [ d2p± p\ + {Mz-M N (l-z)) 2 

" 3 N J (2VT)3 [ p 2 ± + M 2 Z _ M 2 z{l _ z) + M 2 (1 _ z) ] 2 ' 



for both cases N = p or n. The remaining function dj^(z) is then obtained from the relation 
This concludes the evaluation of the fragmentation functions of Figs. [3^i, [3f3. 

Here, we note that another possible channel for a nucleon emission within NJL formalism, 
q — > Dq and subsequent D — > Nq is not included in the current version of the model, as numerically 
the norm for the splitting function of this channel is 2 orders of magnitude smaller than the norm 
of the splitting function for the channel described above, thus proving insignificant. 



III. MONTE CARLO APPROACH TO CALCULATING THE FRAGMENTATION 

FUNCTIONS 

A. Monte Carlo Simulations as an Alternative to the Integral Equation Method 

The Monte Carlo method for describing quark fragmentation process using most notably the 
Lund string model [13] has been long employed to successfully describe the hadronization process 
|44j and has been implemented in very sophisticated event generator frameworks like PYTHIA [45j . 
These frameworks have been refined, extended, and tuned over several decades to cater for the 
needs of the experimental data analysis. Our purpose here is to develop a standalone MC sim- 
ulation software that has the specific purpose of implementing the quark-cascade description of 
the hadronization process with quark splitting functions supplied from an effective quark model, 
particularly NJL-jet model in the current article, and vector meson decays to pseudoscalar mesons. 
This allows for flexibility and simplicity of the software platform, making it readily accessible for 
further development. 

Previously, in the NJL-jet model the fragmentation functions were obtained as solutions of a 
set of coupled integral equations ( [24, [25] ) , 

D h q (z)dz = d h q (z)dz + Y, f d%(y)dy D h Q ( Z -\ ^, (17) 
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where D^(z) denotes the fragmentation function of quark q to hadron h carrying light-cone mo- 
mentum fraction z, and dg(z), d®{z) are the elementary fragmentation functions for the process 
q — > hQ, normalized as Ylh I dq(z)dz = Y1<q J d®(z)dz = 1, thus allowing an interpretation as the 
probability of an elementary process. The sum on the right-hand side is over all possible inter- 
mediate states in the quark cascade ( in our model that includes all the active flavors of quarks 
and scalar antidiquarks). In formulating the integral equations, one assumes that the quark has 
infinite momentum and produces an infinite number of hadrons, which physically corresponds to 
the Bjorken limit. 

We propose to calculate the fragmentation functions using Monte Carlo simulations akin to the 
method described in Ref. [36] using the probabilistic interpretation: D^.(z) dz is the probability 
to emit a hadron h carrying the light-cone momentum fraction z to z + dz of initial quark q in a 
quark-jet picture. The quark goes through a cascade of hadron emissions, where at every emission 
vertex we choose the type of emitted hadron h and its light-cone momentum fraction z (of the 
fragmenting quark) by randomly sampling the corresponding probability densities of the elementary 
fragmentations, dg(z), that are calculated within the NJL model (in general these can be calculated 
in any effective quark model). We keep track of the flavor and the light-cone momentum fraction 
of the initial quark left to the remaining quark, also recording the type and light-cone momentum 
fraction of the initial quark transferred to the emitted hadron in each elementary fragmentation 
process. We stop the fragmentation chain after the quark has emitted a predefined number of 
hadrons, N Links- (Alternatively, one could stop the chain after the remnant quark in the cascade 
has less than a given fraction of the initial quark's light-cone momentum, ZMin-) We repeat the 
calculation Nsims times with the same initial quark flavor q, until we have sufficient statistics for 
the emitted hadrons. We extract the fragmentation functions by calculating the average number of 
hadrons of type h with light-cone momentum fraction z to z + Az, (Ng(z, z + Az)} and expressing 
it as 

/ u \ Y.n Nj}(z,z + Az) 
D h Q {z)Az = (N*(z, z + Az)) = ^ n ^s _g v ( 18 ) 



q\~, yq v->- ■ // - w 

From the construction, it is obvious that the fragmentation functions calculated using the 



integral equations, Eq. (17), should be equivalent to those calculated using the MC method in the 
limit N Links co and Nsims oo. The plots in Fig. [4jshow that the solutions for fragmentation 
functions from both methods are indeed equivalent with a large enough number of emitted hadrons 
within statistical errors. It follows from Eq. (17) that the solution to the integral equations behaves 
as zDq(z) — > const as z — > 0. This behavior originates in the assumption made by Field and 
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Feynman in the original jet model [23J that the total fragmentation function can be expressed 
as an infinite product of elementary fragmentation functions, which allows one to formulate the 
integral equations. The deviations of the MC solution from the solution of the integral equations for 
very small z is because in the MC calculation we always use a finite number of hadrons emissions, 
although in Fig. [4] this number was taken to be large enough so as to demonstrate the equivalence 
to the integral equations. 



0.9 



0.4 



Integral Equations 

X Monte Carlo Simulations 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



1.0 



FIG. 4. Comparison of the solutions for quark fragmentation function (z) in NJL-jet model with only 



nonstrange pseudoscalar mesons calculated from integral equations, Eq. (17) and MC simulation 



MC also allows us to study the dependence of the resulting fragmentation functions on the 
number of hadrons emitted by the quark in the cascade, which could well be relevant to many 
medium-energy experiments. Figure [5] shows that the solution for zD^ + (z) with NLinks = 1 
[equivalent to the elementary fragmentation function zd^ + (z)\ is peaked at z ~ 0.8. As we increase 
the number of emitted hadrons, the solution increases in the low z region because of the hadrons 
emitted further in the quark jet, where the fragmenting quark typically has a small fraction of 
the initial quark's light-cone momentum. We can readily see that the solutions saturate after 
including only a few emitted hadrons, where there is virtually no difference between solutions with 
N Links = 8 and NLinks = 20, and the discrepancy to the solution of the integral equations only 
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occurs at extremely small values of z, approaching the limit zD(z) — > const in the case of large 
number of links. Thus we can reliably use the solutions of MC simulations with NLinks > 8. 




FIG. 5. The dependence of the solutions for zDJ on Nn n k s . 



B. Including the Resonance Decays 



In the current version of the NJL-jet model, we included the production of vector mesons, among 
the other hadrons directly emitted by the quark cascade (so called "primary" hadrons). The vector 
mesons considered have an extremely short lifetime and decay quickly, thus in an experimental 
setup one usually detects only their decay products ( "secondary" hadrons), most often pions 
and kaons. Schematically, the process is depicted in Fig. [6j Consequently, to best describe the 



experimentally measured fragmentation functions of a hadron h from the Eq. (18), we should not 
only consider the number of primary hadrons with a given range of light-cone momentum fraction 
of the original quark, N^(z 7 z + Az), but also add the number of hadrons h satisfying the above 
criteria that originate from the decays of primary vector mesons (or in general other hadronic 
resonances). This is accomplished using the dependence of the decay probability of a hadronic 
resonance h on the fractions of its light-cone momentum, z\ to zi + dzi, z n to z n + dz n ; ^ Zi = 1, 
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carried by the decay products hi, ...h n , denoted as dP 1 " , ' ln (zi, z n ). First we perform the MC 
simulations described in the previous section and record all the produced (primary) hadrons along 
with the fractions of the initial quark's light-cone momentum. We calculate the fragmentation 



functions without the decays using the formula in Eq. (18). Second, we consider each encountered 
resonance particle h with a momentum fraction of the initial quark z and its possible strong decay 
channels, randomly selecting one according to the corresponding relative branching ratio. Then 
we randomly generate the light-cone momentum fractions z\ , . . . , z n carried by the decay products 
hi, ...h n according to the probability dP hl '" hn (zi, ...,z n ). The decaying hadron h is removed from 
the list of the produced hadrons and the decay products hi, h n are added with their respective 
fractions of the initial quark's light-cone momenta zz±, zz n . The fragmentation functions are 



once again calculated using the new list of produced hadrons using the formula in Eq. (18). (In 
general, the decay of a primary resonance can produce another resonance of a lower mass that 
decays itself, so we start the simulation of the decay process from the heaviest resonances present 
and move on to the lighter ones.) 



i 



////// 
t — > — • — > — • — >■ 



Q* Q" 

FIG. 6. Quark- Jet model with the decay of the resonances. 

In the current article, we consider only the strong decays of the vector mesons to two-body 
pseudoscalar meson final states for simplicity. A generalization to include the three or more-body 
final states is straightforward, although it requires sampling the nontrivial phase space factors using 
Monte Carlo techniques. We consider all of the two-body strong decays of <j>, K* , and p mesons 
with the corresponding relative branching ratios as given in the " Review of Particle Physics" |47| . 
For calculation of the two-body decay probability function, let us consider the initial hadron h 
with mass mn, momentum q decaying to hadron hi with a mass rrihi, and a momentum p\ and 
hadron h2 with mass m/j 2 , and momentum p2 ■ We also denote the light-cone momentum fraction 
of the h carried by hi as zi = p7/q~ and the fraction carried by hi as Z2 = P2 /q~ , where trivially 
Zi + Z2 = 1. Thus, the decay probability is a function of only one momentum fraction chosen to 
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be the z\. The dP hl% 2 {z\) is determined as a product of the decay amplitude squared times a 
two-body phase space factor. 

A detailed description of the decay amplitudes and the branching ratios into different channels 
can be calculated using specific models ( for example model Lagrangians for the interaction from 
Ref. |48j). Here, we are only concerned with the dependence of the decay probability on zi, where 
we average over the polarization of the vector meson. Thus, the Lorentz invariance requires that 
amplitude squares depend only on scalar products of the 4-momenta of the particles involved in 
the decay, which in turn are trivially expressed through on-mass-shell conditions in terms of their 
masses. Thus, the only dependence on z\ comes from the two-body phase space factor. Our goal is 
to express the elementary probability for the decay as a function of z\, assuming a constant cj^ lh2 
for the amplitude squared of the decay. For that we integrate over all components of momenta 
pi and P2 with exception of p± , assuming qj_ = and using the light-cone form for the two-body 
phase space factor: 

ihr+hxMf^s _nhiha j„- f _^|*Pl± ^£^P^/ x4r4. 



dP ^ {zi) =cr , dpI j ^^- (M - n - n) (19) 

= |- dz 1 J dl 5 (I - [ziz 2 ml - z 2 m 2 hl - z x m\ 2 \) \ Z2 =i-z 1 (20) 

dzi if ziz 2 m 2 h - z 2 m 2 hl - z x m\ % > 0; z x + z 2 = 1, 
otherwise. 

Here we note that the numerical values for the z\ ranges for various decays exactly match those 
shown in the plots in Fig. 2 of Ref. [19]. We determine the C^ lh2 such that the total probabilities 
j jph^hi of the decays to different pairs {hih 2 } relate as the corresponding relative branching 
ratios and the normalization condition ^{/n/i 2 } / dP h ^ hlh ' 2 = 1 is satisfied. 



IV. RESULTS 

The plots in Fig. [7] depict the elementary fragmentation functions, zd^(z) from a u quark 
for pseudoscalar, vector meson and nucleon emission with the mass of scalar antidiquark Mjj = 
0.687 GeV from Ref. [50] and the following quark-hadron couplings (calculated in Appendix |C|) : 
9-wqQ = 3.15, g Kq Q = 3.387, g pqQ = 1.5, §K*qQ = 1-91, and g^ qQ = 2.29. 

In order to compare our results with experimental measurements or empirical parametrizations 
we evolve them at next-to-leading order (NLO) from our model scale Qq = 0.2 GeV 2 using the 
software from Ref. |51j . The details of determination of the model scale are given in [25J. 



15 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

z 



FIG. 7. Elementary fragmentation functions for u quark, zd%, for all included emission channels. 

The solutions for the favored and unfavored fragmentation functions from an u quark to primary 
hadrons (without any resonance decays), evolved at NLO to a typical scale Q 2 = 4 GeV 2 , are shown 
in the plots in Fig. [8j Here, we note that the fragmentation to nucleons is comparable to the case 
of vector mesons. 




0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 

z z 



(a) (b) 

FIG. 8. The solutions for (a) favored and (b) unfavored fragmentation functions from u quark to various 
hadrons evolved at NLO to Q 2 = 4 GeV 2 . 

The plots in Fig. [9] depict the total light-cone momentum fractions of the initial quark carried by 
the emitted hadrons of different species, calculated at model scale of Qq = 0.2 GeV 2 . (As discussed 
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in Ref. |25j . the momentum and isospin sum rules can only be reliably satisfied at model scale, as 
the NLO QCD evolution kernels have known singularities in the low z region, in practice leaving 
the values of the fragmentations for z < 0.01 undetermined.) Here we compare the fractions for 
elementary fragmentation functions, the solutions without the resonance decays and the solutions 
after the resonance decays. The sum of the fractions from elementary splitting functions calculated 
in the NJL model amounts to about half of the initial quark's light-cone momentum, where the rest 
is carried by the remnant quark. In the quark-jet picture, we sum up the light-cone momentum 
fractions of all hadrons emitted in the chain. This implies that the initial quark transfers all of its 
light-cone momentum to the produced hadrons, which we confirm in our numerical solutions (the 
sum of all the fractions is 1 within numerical errors). This shows that after resonance decays, the 
pions carry approximately 67% and kaons carry about 24% of the initial u (or d) quark's light-cone 
momentum. In the case of an initial s-quark, pions carry approximately 25% and kaons carry 
about 72% of its light-cone momentum. The rest of the momentum is carried by nucleons and 
antinucleons, where the momentum sum rule is satisfied within 0.1%. 



The results for the solutions for fragmentation function zD^ + and zD^ prior to and after the 



vector meson decays are shown in Fig. 10 Here, we see that the inclusion the resonance decay 
slightly changes the shape of the function in mid- to low- z region, bringing it closer to the empirical 
parametrizations of the experimental data from Refs. [12] (HKNS) and [10] (DSS). Similarly, the 



results for zDf~ prior to and after the vector meson decays are shown in Fig. 11 and the results for 



zD^ and zD^ are shown in Fig. 12 For zD^ the agreement with the empirical parametrizations 
from Refs. |12] (HKNS) and [11] (DSS) is reasonable, while our results for zD^ are well below the 
empirical curves. This is because we only include the scalar diquarks in our model, making the 
fragmentation of a u quark to N an unfavored process, while in the empirical parametrizations it 
is usually assumed D^(z) — 2D^ (z) based on naive SU{2) flavor symmetry arguments. In the 
future developments of the NJL-jet model, the inclusion of the axial-vector diquark intermediate 
states in the quark to nucleon splitting process will include a favored channel [depicted in Fig. [2^, 
] for emission of a neutron from a u quark. 
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Here we see the strong discrepancies 



The solutions for the zD^ and zD^ + are shown in Fig. 
in the global fits [12] and [TU] of the experimental data that illustrates the need for the model 
calculations providing additional insight into the quark fragmentation process. 
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FIG. 9. The total fractions of momenta carried by mesons of type m from the jet of (a) u and (b) s quarks 
calculated using the numerical solutions for fragmentation functions at the model scale Q% — 0.2 GeV 2 . Here 
"splittings" denote the momentum fractions calculated using the elementary splitting functions (zd™(z)), 
"NJL-jet" and "with decays" denote the momentum fractions calculated from solutions for the fragmentation 
functions without and with inclusion of pions and kaons originating from decays of vector meson resonances 
(the corresponding columns arranged from left to right for each hadron). 
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FIG. 10. The solutions for fragmentation function zD^(z) evolved at NLO to Q 2 — 4 GeV 2 . The results 
are compared to the empirical parametrizations of the experimental data from Refs. |12 (HKNS) and [TU] 
(DSS). Here "NJL-jet" and "with decays" denote the fragmentation functions calculated without and with 
inclusion of pions originating from decays of vector meson resonances. The shadows show the uncertainties 
of the empirical functions of Ref. [H] . 
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FIG. 11. Same as Fig. 10 for the case zD£(z). 



V. CONCLUSIONS AND OUTLOOK 



In the current article, we added the vector meson, nucleon, and antinucleon emission channels to 
NJL-jet framework for calculating quark fragmentation functions. We also included the two-body 
decays of the vector mesons to pseudoscalars, which allows for easier comparison of the calcu- 
lated fragmentation functions with the experimental measurements or their parametrizations. We 
employed the Monte Carlo method to obtain the fragmentation functions in a quark-cascade de- 
scription. Here we demonstrated that the Monte Carlo approach to calculating the fragmentation 
functions in NJL-jet framework is a powerful and reliable method. We reproduced the fragmenta- 
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tion functions calculated as solutions of the previously employed integral equations, where only the 
light quarks and pions were included. Moreover, we showed that the MC approach allows for the 
flexibility to surpass the model limitations necessary in formulating the integral equations. That 
is, in the future MC studies we can assume an initial quark carrying only a finite momentum and 
thus emitting a finite number of hadrons. We demonstrated that the medium and low z regions 
of the fragmentation functions are greatly affected by the number of the emitted hadrons, thus 
the finiteness of the quark momentum might have a noticeable effect. The future development 
of the NJL-jet model would also allow an access to the transverse momentum distribution of the 
produced hadrons, thus becoming relevant for the analysis of a large variety of semi-inclusive data. 
The MC approach provides a robust and efficient platform for implementing these and other pos- 
sible extensions of the NJL-jet model that would allow for a much more detailed description of the 
physical picture. 
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A further advantage of the MC approach is in reducing the numerical task in solving for the 
fragmentation functions when including many more channels for emitted hadrons. Here, solving 
the integral equations requires inverting larger and larger matrices, while the MC procedure can 
be drastically sped up by trivially parallelizing the task and solving simultaneously on computer 
clusters. 

The results for the fragmentation functions exhibit only slight changes with addition of the 
new hadronic channels. In particular, vector meson-quark couplings are relatively small in our 
model, lowering their relative contribution to the fragmentation process. The resulting quark 
fragmentation functions exhibit a good agreement with the empirical parametrizations in high z 
region, staying reasonably close to them also in mid and low z regions in part due to the slight 
modifications of the shape of the functions in this region brought by the effects of the vector meson 
decays. On the other hand, the plots in Fig. [9] show that pions, kaons, nucleons, and antinucleons 
include the most dominant contributions in fragmentation process, thus the inclusion of higher 
resonance states is unlikely to be important. A notable underestimation of the fragmentation of 
u quark to neutrons, shown on plots in Fig. [l2p , demonstrates the need to include the axial- 
vector diquark for a more realistic description of the fragmentations to nucleons, which will be 
accomplished in the future developments of the model. 

In this work, we followed our previous NJL-jet calculations and used the LB scheme to calculate 
the regularized fragmentation functions. However, as we have pointed out, the inclusion of vector 
mesons in the NJL model favors the PT regularization scheme, which is free of unphysical decay 
thresholds. Here, we have used this scheme only to assess the vector meson-quark couplings, but 
in future work on extensions of the model we will use the PT scheme consistently throughout. 

The future development of the NJL-jet model using the Monte Carlo framework will allow us 
to study the transverse momentum dependence of the quark fragmentation functions as well as 
polarized fragmentation functions. This can be accomplished within the NJL framework, without 
inclusion of any additional parameters. 

Last, the medium modifications of the fragmentation functions may be essential to our under- 
standing of the semi-inclusive processes on nuclear targets. Medium effects have long been studied 
in the NJL model, yielding a successful description of modifications of nucleon properties in nu- 
clei |22[ W2\ I52j . Thus the model should provide a reliable framework for the calculation of 
medium modifications of quark hadronization process. 
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Appendix A: Flavor Factors for Distribution and Splitting Functions 



The flavor factors given in Table I are calculated from the corresponding SU(3) flavor matrices 
(for details see, e.g., Ref. [33]). With the interaction Lagrangian in the vector channel given in 



Section II B the <f> meson is a pure ss state, i.e. an ideal mixture of flavor singlet (Ao) and octet 
(As). 



TABLE I. Flavor Factors C£ 
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Appendix B: Vector Meson Splitting Functions 



In this appendix we present the details for derivation of vector meson splitting functions. We 
start with the expression in the Eq. ([5]): 

™ n = ^ 2 z f d 2 k T 1 Tr[(^ + Af 1 ) 7 +(^ + M 1 )7 M (^-^ + M 2 ) 7 ,] 

q[Z) 2 9m ^2j (27r)3 2p-(l/z-l) (k 2 - M?) 2 



ltw + ~jr ) U-=p-/*, (fc- P ) 2 =Af| i ( B 1 ; ' 
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where we used the following relation to sum over the vector meson polarizations: 



A==Fl,0 ^ P / 



We split the TV in the Eq. (Bl) into two parts: 



T4 = TV 



+ + MO ( - 7 '*(# -f + M 2 ) 7jU + rvr "^ 2 +M2) 



Here 



Trl 



V 



-Tr 



P 



zl-z) 



'0 + Mi) 7 +(^ + Mj-ftf -f + M 2 ) 7m ] 
{z 2 k\ + (M 2 - (1 - ^)Mi) 2 - 2(1 - z)M 1 M 2 ] 



And 



(B2) 



Tr9 d 



(B3) 

(B4) 
(B5) 



TV2^ = TV [(# + Mi) 7 +(# + Mj 



+ M 2 )p] 



P 



{((1 - z)p 2 + z 2 (4 + MiM 2 )) 2 + 4(Mi - M 2 )V} . 



(B6) 
(B7) 



z(l — z)z 2 

In the above derivations, we used the on-mass-shell condition for the emitted vector meson with 
mass m m and fragmented quark in the frame where pt = 0. Thus, we obtain for the elementary 
splitting function 

jm< \ _ C T 2 [ d2 P± 

a q \ z ) — ~^~9mqQ z 



(2tt) 3 

2 (pi + ((1 - z)M 1 - M 2 ) 2 ) + pig- [{pi - z 2 M x M 2 + (1 - z)m 2 m ) 2 + p\z 2 (M l + M 2 ) ; 

(p 2 + z(s - 1)M 2 + zM 2 + (1 - ^m^) 2 



• (B8) 



Appendix C: Quark-Meson Couplings with Proper-Time Regularization 



Here, we derive the vector meson-quark coupling constants in NJL model using the PT regu- 
larization scheme [37 39J to calculate the quark-vector-meson couplings. 



1 



1 



dr T n ~ x e- rX -> 



1 



l/A? 



dr T n ~ L e 



n-l-rX 



(CI) 



X n (n-iy.J (n-l)\J 1/A 2 uv 
where X denotes the denominators of the loop integral after an appropriate momentum shift and 
Wick rotation, Am and Auy are the infrared (IR) and ultraviolet (UV) cut-offs, respectively. It 
was demonstrated in Refs. |37t |3"B] that Air mimics confinement by eliminating the unphysical 
thresholds for the decay of the hadrons into quarks, while in Ref. [39J it has been shown that this is 
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crucial to describe saturation of the nuclear matter binding energy in the NJL model. Here we fix 
M u = 0.3 GeV as in our previous study |25| and Am = 200 MeV as in previous model calculations 
of Ref. |42] and fit Ajjy to reproduce the experimental value of pion decay constant f n = 93 MeV. 
We then obtain the strange constituent quark mass requiring the calculated kaon mass to reproduce 
the experimentally measured value mx = 495 MeV. We will verify that the pseudoscalar meson- 
quark couplings are close to the values calculated in the LB regularization scheme, indicating that 
it is reasonable to use the PT scheme for the calculation of the vector meson-quark couplings, while 
keeping the LB scheme for the regularization of the elementary fragmentation functions. 



1. Pseudoscalar Meson-Quark Coupling 



Here we follow our previous work of Ref. [25] . The quark-meson coupling constant is determined 
from the residue at the pole in the quark-antiquark t-matrix at the mass of the meson under 
consideration. This involves the derivative of the familiar quark-bubble graph 

U(p 2 ) = 2NJ J TrfoSiCAOTS&Cfc-p)], (C2) 



1 _ (BW?)\ (C3) 



9 2 mqQ ~ V dp 2 

Here Tr denotes the Dirac trace and the subscripts on the quark propagators denote quarks of 
different flavor - also indicated by q and Q, where the meson of type m under consideration has 



mass m m and flavor structure m = qQ. Using the proper-time regularization of Eq. (CI) we can 
calculate the above integral as 

TTt ^ zat ■ f d4ft -k 2 + k-p + M l M 2 

UiP } = 8A W (2^ (k 2 -M 2 m -P?-M 2 )) (C4) 



J\k f^ndr /' 1 dx[l + T(i3l2(X)p 2 ) _^ i2(X)p2))] e -.A 12 (^) ) 
Z7T Jl/Afj V T JO 



(C5) 



where 



A 12 {x,p 2 ) = (1 - x)Ml + xMl - x(l - x)p 2 , (C6) 
B 12 (x,p 2 ) = x(l-x)p 2 + M 1 M 2 . (C7) 



For the quark-meson coupling we have: 

h/Kly T Jo 



V~Zq = ^\ X ^ R % j 1 dxx(l-x)[3 + r(B 12 (x,m 2 )-A 12 (x,m 2 ))] e ^ A ^ m2 \ (C8) 
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2. Decay Constant and Fixing the UV cutoff. 

Here, we fix the Ajjv via the decay constant of the pion, which is obtained in the NJL model 
from the one loop contribution of the meson decay: 



ip^fm = N c g mqQ 
iN c g mq Q 



P 



d 4 k 
(2tt) 4 



Tr[^ 5 S 1 (k)fy 5 S 2 (k - p)]\ t 



N c g m qQ 
Att 2 



J dx ((1 - x)Mi + xM 2 ) (Ei 



Ai 2 {x,m. 



A 2 



Ei 



A 12 (x,m 2 n ) 

A 2 

IV UV 



(C9) 
(CIO) 
(Cll) 



Using the previously set values of M u = 0.3 GeV and Air = 0.2 GeV, the fit to the experimental 
value of f n = 0.093 GeV yields for the UV cutoff Aj/y = 0.703 GeV. This allows us to determine 
the strange constituent quark mass, M s , requiring that the kaon mass determined from the quark 
anti-quark t-matrix pole reproduces the experimental value. This gives M s = 0.539 GeV. The value 
of M s in the PT scheme is only marginally larger than the value used in the previous study |25j . 

( l B) 

which was determined using LB regularization: Mr ' = 0.537 GeV. Then we can confirm that 
the values of the pseudoscalar meson-quark coupling constants calculated in both regularization 
schemes are indeed very close to each other 



& = 3-15, = 3.39, 

= 3 - 16 > 9$% = 3.41. 



(C12) 



3. Vector Meson-Quark Couplings 

The vector meson-quark coupling constant, g m qQ, can be calculated in the familiar manner. 
First, we evaluate the quark-bubble graph in vector channel using the PT regularization scheme; 
here only the transverse component contributes to the corresponding t-matrix for our model La- 
grangian with only the four-point quark-quark interactions [33] 



p2 



(f - v -f-) ^£ dx / v ( MiM2 _ (i _ x)m * ~ xM * + 2x[i ~ x)p2 ^ e ' TAi2{x,p2) ■ 

(C13) 
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Thus we finally obtain for the coupling constant: 



r( T ) 

'\p 2 =™ 2 m 

-—^ dxx(l-x) dr{ MiM 2 - (1 - x)M? - xM$ + 2x(l - x)m 2 m + - e ~ TAia(a,>m ™ ) . 

The masses of the vector mesons can be determined within the NJL framework by fixing the 
quark coupling in the Lagrangian using the experimental value for the p meson. Then the masses 
for the K* and 4> mesons are 



-a _ gW , a 2 fci4l 



m£ JI,) = 0.936 GeV, mf JL) = 1.088 GeV, (C15) 

which are reasonably close to the experimental values, given the simplistic Lagrangian used. In 
this work we choose to use the experimental values for all the vector mesons in calculating the 
couplings and the splitting functions in an attempt to best describe the measured fragmentation 
functions. Then the corresponding values of the couplings for the vector mesons considered in this 
article are 

9 P qQ = 1-5, 9K*qQ = 1-91, g^ qQ = 2.29. (C16) 
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